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SIMULATIONS 
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ABSTRACT 

We present results obtained with a three-dimensional, Lya radiative transfer code applied to a fully 
cosmological galaxy formation simulation. The developed Monte Carlo code is capable of treating an 
arbitrary distribution of source Lya emission, neutral hydrogen density, temperature, and peculiar 
velocity of the interstellar medium. We investigate the influence of resonant scattering on the appear- 
ance and properties of young galaxies by applying the code to a simulated "Lyman-break galaxy" at 
rcdshift z = 3.6, and of star formation rate 22 M Q yr _1 and total Lya luminosity 2.0 x 10 43 erg s^ 1 . 
It is found that resonant scattering of Lya radiation can explain that young galaxies are frequently 
observed to be more extended on the sky in Lya than in the optical. Moreover, it is shown that, for 
the system investigated, due to the anisotropic escape of the photons, the observed maximum surface 
brightness can differ by more than an order of magnitude, and the total derived luminosity by a factor 
of ~ 4, depending on the orientation of the system relative to the observer. 

Subject headings: galaxies: formation — line: formation — line: profiles — radiative transfer - 
scattering 



1. INTRODUCTION 

The Lya line is a very important diagnostic in a wide 
range of fields of astrophysics, not the least of which is 
galaxy formation, providing us with extensive informa- 
tion on redshift, dynamics, kinematics, morphology, etc. 
Three distinct physical processes result in Lya source 
emission in the context of galaxies: First, Lya emission 
due to photo-ionization of hydrogen atoms by UV radia- 
tion from nearby, massive stars and subsequent recombi- 
nations may contribute as much as 10% of the to tal lu- 
minosity of the galaxy (|Partridge fe Peebleslll967f ). Sec- 
ond, part of the potential energy gained by gas falling 
into galact ic potential wells i s converted into cooling 
radiation. iFardal et al.l (|2001f ) find that, at high red- 
shifts, most of this radiation is emitted by gas with 
T < 20 000 K, and consequently ~ 50% in Lya alone. 
Finally, the external, metagalactic UV field, penetrating 
some (in case of damped Lya systems) or all (in case of 
Lyman-limit systems) of the outer parts of galactic hy- 
drogen "envelopes" will also produce some Lya radiation 
through case B recombination. Moreover, this UV field 
can also photo-heat non-self- shielded gas, which su bse- 
quently cools, radiating Lya (jFurlanetto et al.ll2005h . 

Over the past years it has become possible to actually 
resolve obscrvationally these young, Lya emitting galax- 
ies. In several cases, the galaxies have been found to be 
significantly more extended on the sky when observed in 
Lya as opposed to optical ban ds (e.g.. lM0ller fc Warrenl 
119981: iFvnbo et al.ll2001l l2003f ). Due to the complexity 
and diversity of the systems, in order to correctly inter- 
pret observations and make predictions about the prop- 
erties of young galaxies, it is desirable to develop realistic 
theoretical and numerical models. 
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For a nu mber of ideal i zed c ases analytical solution are 
obtainable. iHarringtonl (|1973l ) investigated the emergent 
spectrum of resonantly scattered radiation in the case of 
a highly optically thick slab of finite thickness but in- 
finite extensio n, and uniform temperature and density. 
iNeufeldl (|1990f ) extended this solution to include the pos- 
sibility of the photons being destroyed (a s by dust) and 
injecte d with arbitrary initial frequency. iDijkstra et alJ 
(|2006f ) derived a similar analytical expression for spher- 
ical symmetry, allowi ng for isotropic e xpans ion or col- 
lapse of the gas, and lLoeb fc Rvbickil (|1999f l examined 
the spectrum for an isotropically expanding or contract- 
ing medium with no thermal motion. 

However idealized these configurations may seem, they 
provide valuable and at least qualitative insight into the 
characteristics of young galaxies. Moreover, they offer 
direct means of testing numerical methods. 

The Monte Carlo (MC) method has been used for 
solving radiative transfer (RT) problems since the be- 
ginning of the 1960s. Nevertheless, though concep- 
tually simple, the demand for strong computer power 
until quite recently restricted this technique to deal 
with more or less the same idealized configurations that 
had already been dealt with analytically. Thus, the 
majority of previous attempts to model RT in astro- 
physical situations have been based on strongly sim- 
plified configur ations of the physical parameters. Only 
a few authors (ICantalupo et alJ 120051 iTasitsiorml 120061 : 
IVerhamme et al.ll2006l ) considered more general cases. 

With the aim of predicting the appearance and 
properties of Lya emitting galaxies, we developed an 
MC code capable of treating an arbitrary distribution 
of source Lya emission, neutral hydrogen density, 
temperature, and peculiar velocity of the medium and 
subsequently applied it to a simulated Lyman-break 
galaxy (LBG) from a fully cosmological simulation. 
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2. SIMULATIONS 

2.1. The Cosmological Simulation and Monte Carlo 

Code 

The cosmological simulation of the formation and 
evolution of an individual galaxy was performed us- 
ing the N-body/hyd r odyna mical Tre eSPH code of 
iSommer-Larsen et al.l (|2003l see also iSommer-Larsenl 
()2006f) ). The system becomes a Milky Way/M31-like disk 
galaxy at z = and is simulated using ~ 2.2 x 10 6 par- 
ticles in total, comprising only smoothed particle hydro- 
dynamic (SPH) and dark matter (DM) particles at the 
initial redshift Zi = 39. The masses and gravity soften- 
ing lengths of SPH and star particles are 9.9 x 10 4 h~ 1 MQ 
and 200h~ 1 pc, respectively (h = 0.65). For the DM par- 
ticles the corresponding values are 5.7 x 10 5 /i^ 1 M Q and 
370/i _1 pc. The minimum SPH smoothing length in the 
simulation is about l'Sh^ 1 pc. For the purposes of the 
present study, a single simulation output, at z = 3.6, was 
chosen. 

The MC code, used to propagate the photons through 
the medium, by and large rese mbles those recently de- 
veloped by other authors (e.g.. ICantalupo et alj 120051 ; 
lTasitsiomi([2006l ). Since it is grid-based — the number of 
cells typically being 512 3 — the physical parameters of 
interest are first interpolated from the SPH particles to 
the cells of the computational box. These parameters are 
the Lya emissivity, the temperature T, the density riHi 
of neutral hydrogen, and the three-dimensional velocity 
field Vbuik of the gas. 

2.2. Determination of Source Lya Emission 

The emission cell of a given Lya photon is found by 
setting the probability of being emitted from a given cell 
equal to the ratio of the luminosity in that cell to the 
total luminosity L to t- The photon is then injected in the 
line center (in the reference frame of the fluid element) 
from a random point Xj in the cell. The frequency v of 
the photon is parametrized through x = (y — i/q)/ A.i/q, 
where vq = 2.466 x 10 15 Hz is the line center frequency 
and A^d = {vt\Jc)vQ is the Doppler frequency, with 
u t h = {2k bT '/ma) 1 / 2 being the thermal atom velocity 
dispersion (times v2) and the rest of the variables hav- 
ing their usual meanings. In terms of this variable the 
injected photon obviously has a frequency of x = 0. The 
initial direction ft; of the photon follows an isotropic 
probability distribution. 

2.3. Propagation of the Radiation 

The optical depth r covered by the photon before it is 
scattered is drawn randomly from the probability distri- 
bution e~ T , and subsequently converted into a physical 
distance r = r/n-HiOx, where the physical parameters are 
given by the present cell. The cross section a x is given 
by a Voigt profile, i.e. the convolution of the Lorcntzian 
natural line profile and the Gaussian thermal broadening 
of the atoms, resulting in 

<r x = fi2 H(a,x), (1) 

where / 12 = 0.4162 is the Lya oscillator strength, and 
a f + °° e~ y2 

BM -*L («-»)*+*■ » <2) 



is the Voigt function with a = Ai/l/2Az^d the ratio be- 
tween the natural line width Az/l = 9.936 x 10 7 Hz and 
the Doppler width. Since eq. @ is no t analytically inte - 
grable, we use the analytic fit given bv lTasitsiomil (|2006f) . 

If the final position xj of the photon is outside the 
original cell, the photon is placed at the point Xj nt of 
intersection with the face of the cell and the above cal- 
culation is redone with the parameters of the new cell, 
an optical depth r' = r- |x int - Xj| (nm er x ) prev .ceii, and 
a frequency Lorentz transformed to the bulk velocity of 
the new cell. This procedure is repeated until cither the 
originally assigned value of r is spent, and the photon is 
scattered, or it leaves the computational box. 

2.4. Scattering 

Except for a small recoil effect, the scattering is coher- 
ent in the reference frame of the atom. However, to an 
external observer, the non-zero velocity u = v a t om /wth of 
the atom will shift the frequency of the photon. In the 
directions perpendicular to fi^, the velocities u_li.2 will 
follow a Gaussian distribution. When x ~ 0, the photon 
barely diffuses spatially. Only when it has diffused suffi- 
ciently far in frequency space, will it be able to make a 
large journey in real space. To skip these non-important 
core scatterings and thus accelerate the co de, if \x\ is less 
than some critical value £ C rit, following iDiikstra et al.l 
(|2006h uxi,2 is drawn from a truncated Gaussian so as 
to favor fast moving atoms and artificially push the pho- 
ton back in the wing. x cr it is determined according to 
aTQ in the given cell. If aro < 1, a proper Gaussian is 
used; otherwise we find that x cr it = 0.02 e^ ln aT °, where 
(£, x ) = (0.6, 1.2) or (1.4,0.6) for ar < 60 or ar > 60, 
respectively, can be used without affecting the final re- 
sult. The effect of this acceleration scheme is a speed-up 
of several orders of magnitude. 

Due to the resonance nature of the scattering event, the 
velocity u\\ parallel to fjj depends o n x, and is generated 
following [Zheng fc Miralda-Escude! (|2002f) . 

The final frequency x' of the scattered photon (in the 
reference frame of the fluid element) depends on direc- 
tion in which the photon is scattered. For scattering 
in the line center, transitions to the 2Pi/ 2 state results 
in isotropic scattering, while the 2P 3 / 2 transition causes 
some polarization, resulting in a probability distribution 
W{9) cx 1 + | cos 2 9, where 9 is the angle between ft, 
and the outgoing direction ft/ ((Hamilton! |1940T) . Since 
the spin multiplicity is 2 J + 1, the probability of be- 
ing excited to the 2P 3 / 2 state is twice as large as be- 
ing excited to the 2P 1 / 2 state 3 . For scatterings in the 
wing, polarization for tt/2 scattering is maximal, result- 
ing in a dipole distribution W{9) oc 1 + cos 2 9 (jStenflol 
Il980f ). The transition from core to wing scattering is 
given by the value of x where the Lorentzian starts dom- 
inating the Gaussian profile, i.e. where a/irx 2 ~ y/ire~ x , 
or x ~ 3 (the exact value is not crucial). In all cases, 
the scattering is isotropic in the azimuthal angle <fi. In 
the observers frame, the final frequency is then given by 
x' = x — uii + hf ■ u + g(l — Hi ■ ft/), where the factor 
g = hpii/Q I mucvth (|Fieldl ll~959( ) accounts for the recoil 
effect. 

3 For the environments produced here, transitions to the 25 state 
and subsequent destruction of the photon can be neglected. 
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2.5. Observable Surface Brightness Maps 

Following the above scheme, the photons are propa- 
gated through the medium, one by one, until they escape 
the computational box. For each scattering (as well as 
the emission) of a photon, the probability that the pho- 
ton will escape in the direction of six virtual observers 
situated in the negative and positive directions of the 
three principal axes is calculated. This probability is 
added as a weight to a three dimensional array (the two 
spatial dimensions of the projected image of the galaxy, 
plus a spectral dimension for each pixel). The array is 
finally collapsed along the spectral dimension and along 
the spatial dimensions to yield the image and spectrum, 
respectively, that an external observer would see. The 
contribution of each pixel element to the surface bright- 
ness (SB) is then 



loq(SB/erq s"' cm" 1 ) 



2.7 -1.0 0.8 



SB 



Li 



1 1 



plx ff-i J2 n 



E 



'W(0), (3) 



where n p h is the number of photons, is the luminosity 
distance to the observer, f2 p j x is the solid angle subtended 
by the pixel, and r csc is the optical depth of the gas lying 
between the scattering event and the edge of the compu- 
tational box in the direction of the observer. 

The MC code was tested on various simple configura- 
tions for which analytic solutions exist, discussed in the 
introduction. Our code exquisitely passes all tests. The 
results of these tests will be presented in a future paper. 

3. RESULTS 

The MC code was applied to a proto-galaxy at z = 3.6, 
consisting of two small, star-forming "disks" separated 
by approximately 2 kpc, on one of which the computa- 
tional box is centered, and a third more extended disk 
of lower star formation rate (SFR), located about 15 kpc 
from the center. The star-forming regions arc embedded 
in a significant amount of more diffuse, non-star-forming 
Hi gas in a 10-15 kpc thick, sheet-like structure, taken 
to constitute the x-y plane. The total SFR of the sys- 
tem is 22 M Q yr" 1 . Observ ed LBGs have SFRs in the 
range 10-1000 M yr" 1 fe.g.. lRigopoulou et al.|[2006h . so 
the simulated galaxy corresponds to a small LBG. The 
Lya emissivity is produced according the three differ- 
ent scenarios described in the introduction. Specifically, 
the luminosity originating from the Hn in the vicinity 
of massive stars (accounting for approximately 90% of 
Ltot) is determined fo llowing iFardal et all (|2001f ). using 
the code Starburst99 (|Leitherer et al.lll999f ) to yield the 
Lyman continuum and assuming a Miller-Scalo initial 
mass function, a mean Lyman continuum photon energy 
of 1.4 Rydberg, and that 0.68 Lya photons are emitted 
per photonionization. The length of the box used for the 
MC calculations is 50 kpc (physical), and <1l ~ 34 Gpc. 

Fig. [1] shows the results obtained (all with a resolu- 
tion of 512 3 cells), viewed from two different directions 
— from the negative x- and z-direction, corresponding to 
an edge-on and a face-on view of the sheet-like structure, 
respectively. Upper panels assume that the gas is opti- 
cally thin to the Lya line, lower show the corresponding 
results with resonant scattering included. 

The effect of the scattering is incontestable: although 
the original constellation of the principal emitters is still 




Fig. 1. — Bolometric surface brightness map of a simulated 
galaxy at z = 3.6. Left and right column shows the system when 
viewed edge-on and face-on, respectively. The top panel displays 
the galaxy as if the Lyo radiation was able to escape directly, with- 
out scattering. The bottom panel shows the effect of the scattering. 



visible, the surface brightness distribution is clearly much 
more extended. Moreover, we notice the effect of the 
viewing angle. Qualitatively, we expect the photons to 
escape more easily from the face of the sheet than from 
the edge and hence that the system should have a higher 
surface brightness than when viewed edge-on. Here, this 
anticipation is quantified: Fig. [2] shows the azimuthally 
averaged SB profiles. To allow for direct comparison with 
observations, the profiles are also shown not including 
the luminosity of the emitter lying in the outskirts of 
the image and smoothed with a PSF corresponding to 
a seeing of 0.8". The maximum SB of the x-y plane is 
found to be 6.1 x 10 -2 erg s _1 cm -2 , or ~ 14 times higher 
than that of the y-z plane. The average SB, from which 
the total luminosity is usually derived assuming isotropic 
emission, is 3.8 times higher. 

Finally, we show the emergent spectrum. Though 
not as clear for the edge-on view, both profiles ex- 
hibit the characteristi c double-peaked profile (see, e.g., 
IVenemans et al.lfeOOSf ). Obviously, this is the result of 
the high opacity for photons near the line center; diffus- 
ing to either side of vq quickly decreases r so as to make 
escape more probable. Furthermore, both profiles imply 
a net inward velocity of the gas; for infalling gas, red 
photons are shifted into resonance, while blue photons 
escape even more easily, thus enhancing the blue peak 
and diminishing the red. 

Fitting a "Neufeld profile" to the observed spectra can 
give us an idea of the intrinsic characteristica of the sys- 
tem. Unfortunately, there is a degeneracy in that the 
profile is dependent only on the parameter oro, Assum- 
ing some temperature, e.g. 10 4 K, representative of most 
of the Lya emitting gas, one could in principle deduce 
the equivalent column density 2Vhi- Alternatively, if iVni 
is obtainable due to, e.g., the presence of a background 
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Obs. R-band 
Obs. Lya 
Edge — on, true SB 
Edge — on, 0.8' 
Face-on, true SB 
Face-on, 0.8" 





Fig. 2. — Left: Bolometric surface brightness (SB) profiles of the galaxy when viewed edge-on (orange) and face-on (green). Both the 
true (dotted) and the smoothed (solid, see text) profiles are displayed. Also shown are the SB profiles of the galaxy LEG02138-29 (SFR 
~ 15 Mq yr — 1 if extinction is negligible, [Fvnbo ct al. 2003) in Lya (blue) and in the il-band (red, normalized to the maximum observed 
Lya SB). In particular the SB of the y-z plane nicely reproduces the observed SB. Left axes measure the SB at the source, while right axes 
measure the SB observed at z = 0. Right: Spectrum of the emergent radiation, clearly displaying the characteristic double-peaked profile. 
The enhanced blue peak indicates a net inward bulk velocity of the gas. 

quasar, constraints can be put on the temperature. 

A bulk rotational motion of the gas will also alter the 
profile. Since in fact a full spectrum is obtained for every 
pixel element, this effect can be studied through long-slit 
spectroscopy. 

4. DISCUSSION AND CONCLUSION 

The present MC calculations do not include the effect 
of dust. Effectively, dust will act as a photon sink and an 
extra scattering possibility. Since, on average, each pho- 
ton scatters ~ 4 x 10 8 times and travel a distance of ~ 40 
kpc before escaping (determined from a non-accelerated 
run of the code with a few 10 3 photons), approximately 
half of which is in the high density (tihi ^ 0.1 cm -3 ), 
cooler (T ~ 10 4 K) regions, even a small amount of dust 
may be expected to be capable of causing a significant 
decrease in the observed intensity. However, it is not 
clear to which extend the dust will affect the observa- 
tions. If the medium is clumpy, the ratio of Lya to con- 
tinuum radiation m ay in fact be increased dNeufeldll 1 99ll : 
lHansen fe Ohll2006f ). Moreover. lTasitsiomil (|2006h argued 
that dust acting as catalysator for hydrogen molecules 
will lower tihi, making the medium more transparent. 
We will study these effects in future work, implement- 
ing a realistic model of the dust based on the 8 different 
metal species kept track of by the cosmological simula- 
tion. 



As a very recent improvement of the cosmological 
simulation, a post-processed RT scheme of UV ra- 
diation from star-forming r e gions was developed by 
iRazoumov fe Sommer-Larsenl (|2006l ). They found that 
up to 5-10% of the ionizing photons escape these re- 
gions at z ~ 3.6. However, a very preliminar analysis 
indicates that the results presented in this work are not 
significantly changed by the inclusion of H-ionizing pho- 
ton RT. The quantative effects of this on Lya RT will 
also be discussed in a forthcoming paper. 

The developed Monte Carlo code has reproduced qual- 
itatively and quantitatively the observation that young 
galaxies often appear significantly more extended on the 
sky in Lya than in the optical. Furthermore, we inves- 
tigated the impact of the viewing angle on the observed 
surface brightness. Future simulations of a large statisti- 
cal sample of galaxies, taking properly into account dust 
and H-ionizing UV photon radiative transfer, will allow 
us to learn more about the enigmatic Lya emitters. 

The authors are very grateful to J. Fynbo and 
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for quick response to questions, and to J. Schaye for use- 
ful comments. We also thank the anonymous referee for 
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Dark Cosmology Centre is funded by the DNRF. 
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